% -----------
% defining the variables and parameters
% -----------
close all;
var y pai r rr z g epsr ym paim rm;

varexo etag etaz etar;

parameters
beta paistar rstar tau kappa phii phipai phiy old_ylead old_ylag old_rpai
new_ylead new_ylag new_rpai
rhoz rhog rhor
lambda elas mu teta corinov gammap gammac;

 
% -----------
% setting the parameter values
% -----------


%Slope IS-Curve
tau=1/2;

%Slope Phillipscurve
kappa=0.5;

%Steady State Values
lambda=0.8;
elas=3;
teta=0.75;
rstar=2;
beta=(1+rstar/100)^(-1/4);
paistar=4;
mu=0.2;
gammap=0.5;
gammac=0.5;
%Policy Parameters
phii=0;
phipai=0.8;
phiy=0;

%AR(1)
rhog=0.7;
rhoz=0.7;
rhor=0.7;

%Standard errors

sdz=1;
sdg=0.38;
sdr=0.31;
corinov=0.3;

model(linear);

y=((1-(lambda*elas/((1-lambda)*(1+mu)))*(1+gammac*mu/(1+elas*(1-gammac))))
/(1+gammac-(lambda*elas/((1-lambda)*(1+mu)))*(1+(gammac*(1+2*mu))/(1+elas*(1-gammac)))))*y(+1)
+((gammac*(1-lambda*elas/((1-lambda)*(1+elas*(1-gammac)))))
/(1+gammac-(lambda*elas/((1-lambda)*(1+mu)))*(1+(gammac*(1+2*mu))/(1+elas*(1-gammac)))))*y(-1)
-((1-gammac)/(1+gammac-(lambda*elas/((1-lambda)*(1+mu)))*(1+(gammac*(1+2*mu))/(1+elas*(1-gammac)))))*(r-pai(+1))+g;

pai=(gammap/(1+(1+rstar/100)^(-1/4)*gammap))*pai(-1)
+((1+rstar/100)^(-1/4)/(1+(1+rstar/100)^(-1/4)*gammap))*pai(+1)+
((1-teta)*(1-(1+rstar/100)^(-1/4)*teta)/teta)*((1/(1-gammac))+(elas/(1+mu)))*(y-z)
-((1-teta)*(1-(1+rstar/100)^(-1/4)*teta)/teta)*(gammac/(1-gammac))*(y(-1));


r=phii*r(-1)+(1-phii)*(phipai*pai(+1)+phiy*y)+epsr;

z=rhoz*z(-1)+etaz;

g=rhog*g(-1)+etag;

epsr=etar;

ym=y;

rr=r-pai;

paim=paistar+4*pai;

rm=rstar+paistar+4*r;

end;
  





steady;




%shocks;
%var etaz;
%periods 1;
%values 1;
%end;
%simul(periods=50);
%rplot y; 
%rplot pai;
%rplot r;





%-----------
% estimating the model
% -----------

shocks;
var etag; stderr 0.38;
var etaz; stderr 1.00;
var etar; stderr 0.31;
var etag,etaz = (0.3*0.05*0.035);
end; 

%gammac,           0.500,  0.000, 1.000,  beta_pdf,       0.200,  0.100;
%elas,             2.000,  0.100, 5.000,  gamma_pdf,      3.000,  0.500;
%elas,             2.000,  0.100, 5.000,  gamma_pdf,      2.000,  0.500;
%
%Priors are baseline estimates of CGG"
%gammap,           0.500,  0.000, 1.000,  BETA_PDF,         0.5, 0.28867513459481;
%lambda,           0.500,  0.000, 1.000,  beta_pdf,       0.500,  0.100;
%elas,             3.000,  0.500, 5.000,  gamma_pdf,      3.000,  0.500;
%gammap,           0.500,  0.000, 1.000,  BETA_PDF,       0.500,  0.200;
%gammac,           0.500,  0.000, 1.000,  BETA_pdf,       0.500,  0.200;
%gammap,           0.500,  0.000, 1.000,  BETA_PDF,       0.500,  0.200;



estimated_params;
stderr etag,      0.380,  0.100, 0.600,  inv_gamma_pdf, 0.380,  0.200;
stderr etaz,      1.000,  0.420, 3.570,  inv_gamma_pdf, 1.000,  0.520;
stderr etar,      0.310,  0.130, 0.800,  inv_gamma_pdf, 0.310,  0.160;
corr etaz,etag,   0.000, -1.000, 1.000,  normal_pdf,     0.000,  0.400;
lambda,           0.500,  0.000, 1.000,  beta_pdf,       0.350,  0.100;
gammac,           0.500,  0.000, 1.000,  BETA_pdf,       0.500,  0.200;
elas,             3.000,  0.500, 5.000,  gamma_pdf,      3.000,  0.500;
phipai,           0.500,  0.001, 3.850,  gamma_pdf,      0.500,  0.200;
phii,             0.200,  0.000, 0.999,  beta_pdf,       0.500,  0.200;
phiy,             0.150,  0.001, 2.000,  gamma_pdf,      0.250,  0.150;
paistar,          4.000,  0.900, 6.910,  gamma_pdf,      4.000,  2.000;
rstar,            2.000,  0.490, 3.470,  gamma_pdf,      2.000,  1.000;
rhog,             0.700,  0.000, 1.000,  beta_pdf,       0.700,  0.100;
rhoz,             0.700,  0.000, 1.000,  beta_pdf,       0.700,  0.100;

end;


								
										
varobs ym paim rm;																		

%estimation(datafile=data,lik_init=1,mode_compute=0,mode_file=prevolcker_4_mode,load_mh_file,mode_check,bayesian_irf,irf=30,nograph,presample=0,mh_replic=0,mh_jscale=0.3,mh_nblocks=3); 	

estimation(datafile=data,lik_init=1,mode_compute=1,mode_check,bayesian_irf,irf=30,nograph,mh_replic=150000,mh_jscale=0.3,mh_nblocks=3); 	
%estimation(datafile=data,lik_init=1,mode_compute=0,mode_file=prevolcker_4_mode,load_mh_file,mode_check,bayesian_irf,irf=30,nograph,presample=0,mh_replic=200000,mh_jscale=0.3,mh_nblocks=3); 	

(1-gammac)/((1-lambda/(1-lambda)*elas/(1+mu)*(1+gammac*mu/(1+elas*(1-gammac))))+(gammac*(1-lambda/(1-lambda)*elas/(1+elas*(1-gammac)))))


periods=100000;
stoch_simul(order=1,irf=1,simul);

